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O i Abstract. As a generic model for liquid- vapor type transitions in random 

, ' porous media, the Asakura-Oosawa model for colloid-polymer mixtures is studied 

in a matrix of quenched spheres using extensive Monte Carlo (MC) simulations. 

e Since such systems at criticality, as well as in the two-phase region, exhibit lack 

, of self-averaging, the analysis of MC data via finite size scaling requires special 

I . care. After presenting the necessary theoretical background and the resulting 

'"^ ' subtleties of finite size scaling in random-field Ising-type systems, we present data 

P^ ' on the order parameter distribution (and its moments) as a function of colloid and 

Q ' polymer fugacities for a broad range of system sizes, and for many (thousands) 

QJ , realizations of the porous medium. Special attention is paid to the "connected" 

and "disconnected" susceptibilities, and their respective critical behavior. We show 
that both susceptibilities diverge at the critical point, and we demonstrate that 
this is compatible with the predicted scenario of random-field Ising universality. 
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1. Introduction 

Understanding the behavior of fluids that undergo a liquid- vapor phase transition in 
the bulk (or, equivalently, of binary mixtures undergoing bulk phase separation), is 
still rudimentary when one considers the confinement of such systems in mesoporous 
materials, such as porous glasses or silica gels [1]. Such amorphous materials form a 
highly irregular, interconnected, three-dimensionally percolating network, and liquid- 
vapor type transitions of fluids conflned between the walls of such networks are widely 
observed and of practical importance [2]. However, the precise nature of the Hquid- 
vapor critical point of such systems is still only partly understood [3-14]. While 
de Gennes [7] has presented a simple argument that the critical behavior of fluids 
in such random media can be mapped onto the random-fleld Ising model (RFIM) 
[15-17], this prediction could be conflrmed neither by experiments [3-6], nor by 
numerical calculations (MC simulations [8,9,11], or density functional theories [12-14], 
respectively). 

On the other hand, these studies could not point out any flaw in this more than 
twenty year old argument of de Gennes [7] either. In short, the argument of de Gennes 
starts from the well-known phenomenon of capillary condensation [1,18,19]. In an 
inflnitely long slit-pore, the liquid-vapor transition is "shifted" relative to the bulk, 
due to the attractive forces between the fluid particles and the walls. The magnitude 
of the shift depends on the nature of the fluid and the type of walls. In addition, 
if "drying" rather than "wetting" would occur for very thick slits, also the opposite 
effect of "capillary evaporation" may take place [1,18,19]. In general, the chemical 
potential /icoox(£') where liquid and vapor coexist inside a slit pore, differs from 
the bulk coexistence chemical potential ^cocx(oo), and this difference depends on the 
width D of the slit pore. In an irregular interconnected pore network, the local pore 
diameter at position f fluctuates randomly around some average value. As a result, the 
local chemical potential ^(r) where phase coexistence would occur, will also exhibit 
(quenched) fluctuations around some average value (the fluctuations are quenched 
because the structure of the porous network does not change over time) . The analogy 
to the RFIM is readily seen when the fluid is described as a lattice gas, since the latter 
is isomorphic to the Ising ferromagnet. In terms of the Ising ferromagnet, quenched 
random fluctuations in the local chemical potential, become isomorphic to quenched 
external fleld variables hi, with hi a random variable acting on the spin at the z-th 
lattice site (which is precisely the random-fleld Ising model). This reasoning is also 
easily carried over to binary fluid mixtures [7]. 

In the present work, we contribute to the clariflcation of this problem, by 
presenting extensive MC data for a particularly simple model, namely the Asakura- 
Oosawa (AO) model [20, 21] of colloid-polymer mixtures, inside a quenched random 
porous medium. The AO model is known to capture bulk experimental observations 
very well (by bulk we mean in the absence of any porous medium), including phase 
separation [22,23] and interfacial properties [24]. Computer simulations [25-27] have 
shown that bulk phase separation in the AO model, which occurs for sufficiently large 
polymers at sufficiently high polymer fugacity, belongs to the universality class of the 
Ising model [28-31]. In addition, the standard predictions for capillary condensation 
in slit pores [18,19,32] have been well conflrmed for this model [33-36]. Following our 
previous work [37,38], we now consider the more complex problem of the AO model 
inside a random porous medium. The porous medium is obtained using an "easy" 
recipe: we simply distribute a set of obstacles (spheres) at random positions in the 
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simulation box. Once the spheres have been positioned, they remain "fixed", to mimic 
the quenched nature of the medium. Next, the AO model is inserted into the medium, 
and its phase behavior is studied. In particular, we will focus on (appropriately 
constructed) "susceptibilities" of the form [(■)^] — [(■)]^, with (•) the conventional Gibbs- 
Boltzmann thermal average, and [■] an average over many different realizations of 
the quenched obstacles. Of course, in the absence of the porous medium, any such 
"susceptibility" is trivially zero. However, in its presence, the analogy to the random- 
field Ising model implies that such quantities will actually diverge at the critical point, 
and will do so with a characteristic critical exponent 7. 

The outline of our paper is as follows. In Section[2l we recall in detail the necessary 
background of finite size scaling in the Ising and the random-field Ising models, and 
we discuss how these techniques may be carried-over to fiuids with quenched disorder. 
Next, in Section [Sj we define the AO model, explain how this model may be extended 
to also capture quenched disorder, and we describe our simulation method. The results 
are presented in Section IH and we end with a discussion, conclusion, and summary in 
the last section. 

2. Finite size scaling in the Ising model, the random-field Ising model, 
and related models 

2.1. Ising model 

We first consider the pure Ising model, i.e. without any random field, and discuss 
how finite size scaling can be used to extract the critical properties of this model. 
To be specific, we consider a (nearest-neighbor) Ising ferromagnet on a hypercubic 
d-dimensional lattice of linear dimension L and periodic boundary conditions 

Hlsing = -j'^SiSj - H^Si, Sj = ±1, (1) 

with J the exchange constant and H an uniform external magnetic field. Defining the 
instantaneous magnetization per spin s as 

*^;^I]**' (2) 

i 

the object of interest is essentially the distribution 

Pl{s)=Pl{s\T,H), (3) 

defined as the probability to observe a magnetization per spin s, in a system of size 
L, at temperature T and field strength H. Basic observables of interest follow from 
the moments (s'^) = J_^ s'^PL{s)ds of the distribution, where (•) is a conventional 
Gibbs-Boltzmann thermal average. For instance, the average magnetization per spin 
can be written as 

m(T,F) = (,s), (4fl) 

and for the susceptibility x we obtain 

x{T,H) = ^^L''{{.')-{sf), (46) 

where the factor ksT has been absorbed in the definition of x, with fee the Boltzmann 
constant. We emphasize that these expressions must be used with care when Pl(s) 
is bimodal. This happens, for example, at low temperature and H = 0, since then a 
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spontaneous magnetization exists, which may be positive or negative. Consequently, 
Pl{s) has two peaks, one at positive and one at negative values. However, blindly 
applying Eq. l(4a|) . one finds that m{T,0) — irrespective of T, which is not really 
desirable. Since the Ising model has spin reversal symmetry, we have Pl{—s) = 
Pl{+s), and so an easy fix is to introduce 

m'iT,H) = {\s\), (5a) 

x'{T,H)=L^'{{s')-{\s\r), (56) 

which are to replace m and x in these cases. The absolute value has the same effect 
as using a modified distribution {Pl{s) + Pl{—s)) /2 with the integration domain 
restricted to s > 0. Clearly, such a modification is reasonable when Pl{—s) = Pl{+s) 
somewhat holds. For very asymmetric distributions, a safer approach is to define 
m and x in terms of peak positions and widths, respectively. This approach was 
successfully applied to the AO model in [27] and will also be used in this work later 
on. Of course, in the thermodynamic limit, all definitions become equivalent, see 
discussion in [39]. 

As is well known, for d > 2, the Ising model has a second order phase transition 
from the (disordered) high-temperature paramagnetic phase, to the (ordered) low- 
temperature ferromagnetic phase, at some critical temperature Tc. In the vicinity of 
Tc, we expect power law singularities [28] 

?Ti(T, 0) oc {—ty (order parameter), (6a) 

x(T,0)«|t|-^ (66) 

^cx|ir^ (6c) 

with t = T/Tc — 1 the reduced distance from the critical point, and ^ the correlation 
length of the magnetization fiuctuations. In the above, /?, 7, and v are critical 
exponents, which characterize the universality class. 

Of course, the divergence of the correlation length cannot be captured in a finite 
simulation box of size L, and so the above power laws are never observed directly. 
The state-of-the-art is to perform several simulations, using a range of system sizes 
L, and to extrapolate the simulation data to L ^ 00 via finite size scaling. In its 
simplest form, finite size scahng is just the statement that, in a finite system at the 
critical point ^ oc L [28]. Ehminating t from Eqs. (|6aj) and l(6c|) . and using ^ oc L, one 
immediately derives the L-dependence of the magnetization order parameter at the 
critical point 

rriL oc L-^/", (7a) 

Similarly, for the susceptibility, one obtains 

XL^L''/'^. (76) 

These equations simply state that, if one performs a simulation at the critical point 
over a range of system sizes, the magnetization should vanish oc L~^''^, and the 
susceptibility should increase oc L''''^ . 

The above scaling laws are quite general, and should hold near any critical point 
where the correlation length diverges as a power law, i.e. conform Eq. l(6c|) . If also the 
hyperscahng relation is obeyed 

7 + 2/3 = i/d, (8) 
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{a)T<T^ A^/.UJ (b)r>7 




Figure 1. Schematic representations of the distribution Pl(,s) below the critical 
temperature Tc (a) and above (b). Since the susceptibility for T 7^ Tc is finite, 
the peak widths w vanish with increasing system size oc L"'*''^, leading to a 
distribution featuring two 5-peaks when T < Tc, and a single (5-peak when T > Tc- 



with d the spatial dimension, it follows that the entire distribution -Pl(s) at H ^ 
scales with L as [40] 

PLis)\H^o^L^^''p(L/(,sLP/n, (9) 

with Pl{s) the magnetization distribution of Eq.(l3]). Here, p{x,x') is a universal 
scaling function, which essentially depends on the universality class, and the scaling 
should hold in the limits ^ -^ 00, L -^ cxd, with L/£^ finite. Using Eq.(l9]), one readily 
obtains the moments 



{\s\)= \s\PL{s)ds = L-^/''k{L/0. 



{3")= s''PL{s)ds = L-'^P'''h{L/0 ik>0), 



(10a) 
(106) 



which also define the scaling functions /. For k = 1, one recovers the average 
magnetization of Eq. l|4ap . and the expected scaling law Eq. l|7ap is correctly reproduced. 
For the susceptibility, however, these moments imply x ^ i^d—211/v ^ consistent with 
Eq. (|7&l) only when hyperscaling holds. For the Ising model, hyperscaling indeed 
holds [28], and so the use of Eq.|[9]) is justified here. 

Hyperscaling also implies a remarkable property concerning the shape of Pl{s) 
at criticality. To see this, note first that, below Tc, there exists a spontaneous 
magnetization. The magnetization may be positive (+) or negative (— ), and so Pl{s) 
features two peaks centered around s = ±m, see Fig. [Ufa). Each of the peaks may be 
approximated by a Gaussian [40] 

P±(s) « L'"\2TrTxr"^ e^V [-(.s T m)2LV(2Tx)] , (11) 

leading to a squared peak width up = (s^) — (s)^ ~ xT/ L'^- In the thermodynamic 
limit, the peaks remain at their respective positions s ±m. At the same time, since 
the susceptibility away from Tc is finite, the peaks also become increasingly narrow, 
eventually converging to a distribution consisting of two (5-peaks. The behavior 
in the disordered region T > Tc follows analogously. In this case, a spontaneous 
magnetization is absent, and Pl{s) is just a single Gaussian, centered around s = 0, see 
Fig.[TlJb). For T > Tc, Phis) thus converges to a single (5-peak in the thermodynamic 
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limit. Mathematically, the shape of the distribution in these two limiting cases can be 
expressed using the cumulant Ui = (s^)/(|s|)^ IJ. After some algebra, one finds that 

lim Ui — I (two (5-peaks), (12a) 

L^oo,T<Tc 

lim Ui=tt/2 (one (5-peak). (126) 

L — *oo, T>Tc 

Precisely at Tc, the behavior of the cumulant is more subtle [40]. In a finite system 
at Tc, the magnetization vanishes oc L~^''^ , see Eq. fZaj) . Hence, Pl{s) still exhibits 
two peaks, at positions oc ±L~^/^ . For the root-mean-square peak width, we obtain 
w'^ = (s^) — (s)^ oc L'^/"'^'^, where now Eq. l(7&|) was used. Comparing the distance 
between the peaks to their widths, we find 

A ^ , P^"'^ "^^f oc L-, (13) 

peak-to-peak distance 

with u) — [jjv — cf)/2 -f jijv. By virtue of hyperscahng one has a; = 0, implying that 
the relative peak width A does not vanish in the thermodynamic limit. Consequently, 
Pl[s) at criticality does not become a superposition of two (5-functions, but instead 
converges to a distribution of two overlapping peaks. The cumulant U\ at T^ differs 
therefore from the off-critical values [40]. By using Eqs. l|10aP and l|10&p . C/* can be 
expressed in terms of the scaling functions as t/* = /2(i/^)//o (i/<^), which is a 
universal function of L/^, and tends to a universal finite constant. In simulations, this 
result is useful since plots of \J\ versus T for various system sizes L will show a common 
intersection point, yielding an estimate of both t/* and Tc (cumulant intersection 
method [40]). 

Z.Z. random-field Ising model: finite size scaling 

The analysis of MC simulation [39, 41] data for systems belonging to the universality 
class of the random-field Ising model (RFIM) [17] has certain subtleties [42,43], when 
one tries to apply finite size scaling methods [28,40]. The source of the problem 
is that the standard hyperscahng relation [28] between critical exponents, which is 
required by Eq.^ [40], does not hold for the RFIM [17,44-47]. Since the presence 
of the random field breaks the spin reversal symmetry, it is necessary to consider also 
the "disconnected" susceptibility Xdis [45], in addition to the standard "connected" 
susceptibility x- 

The RFIM Hamiltonian reads as 

Hrfim = -'^y^ s^Sj ~ ny^Sj- y^/i^s,, Si = ±1, (14a) 

with J and H defined as before. In addition, at each lattice site i, there acts a quenched 
random field hi, which we take to be completely uncorrelated between neighboring 
sites, and with an average of zero 

hi = ±h, [h,] = 0. (146) 

The amplitude h of the random field should be small but finite {h/J <^ 1), but we are 
not concerned with the crossover to the pure Ising model here, and hence disregard 

I Of course, the cumulant is to be calculated for the full distribution. In particular, for T < Tc, one 
should write Pl{s) = {Pj~ + P^)/2. 
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the limit ft, — > 0. Basic observables are again the average magnetization per spin m, 
the connected susceptibility x, and the disconnected susceptibility Xdis 

m{T,H) = [{s)], (15a) 

x{T,H)=L%s')-{sn (156) 

Xdi.(T,i/) = L'^[(s)2]. (15c) 

For the same reason as before, we also introduce 

m'iT,H) = m)], x'iT,H) = L'[{s') - {\s\n (16) 

For the RFIM model, one has to perform the standard Gibbs-Boltzmann thermal 
average (•) for one realization of the random field, followed by an average over M 
different random field configurations [•], whereby M should be large. Note that Xdis is 
simply the fluctuation of the average magnetization (s) between different reaHzations 
of the random fleld. Due to random variations in these flelds, (s) will sometimes be 
negative, and sometimes be positive. In the limit M — > oo, one has [(s)] = 0, of course, 
but the fluctuation [(s)^] — [(s)]^ will generally not be zero, which is essentially what 
Xdis corresponds to. We shall also be interested in the distributions 

PlAs)=PlAs\T,H^), (^ = 1,...,M), (17) 

deflned as the probability to observe a magnetization per spin s, in a system of size 
L, at temperature T and external fleld Hi, for the i-th random fleld realization. Note 
that we allow Hi to vary between different random fleld realizations. Ideally, one 
would like to have M — > oo, but since resources are limited, simulations always deal 
with flnite M. 

Assuming that the RFIM, for small enough h, has a second order phase transition 
at Tc, we expect power law singularities for m, x, and £, as before, but with different 
critical exponents characteristic of the RFIM universality class [17,45]. In addition, a 
power law is expected for the disconnected susceptibility 

Xdis (X \t\-^, (18) 

with a new critical exponent 7 [17,45]. It has been proved rigorously that the RFIM in 
d = 3 dimensions, at low enough temperature, indeed exhibits a nonzero spontaneous 
magnetization [16]. It has not, however, been proved that the second order transition 
assumed above actually exists (also weak flrst order transitions [48], or spin-glass 
type phases [49] have been suggested). Recent MC simulations, however, favor a 
second order transition, albeit that the critical exponents are still not known very 
accurately [50,51]. 

While for the pure Ising model we have the standard hyperscaling relation between 
critical exponents [28], for the RFIM, rather a different relation has been proposed [44] 

7 + 2/3 = i/(d-6i). (19) 

Here, 6 is an exponent which measures the deviation from the standard hyperscahng 
relation; when it is zero, standard hyperscaling is again recovered. Using the further 
result that 9 = 7/1/ [45], it follows that (7 + P)/v = d/2. 

We now discuss flnite size scaling in the RFIM, following Eichhorn and Binder 
[42,43]. Note flrst that the "derivation" of Eqs.lllal) and ||l6l still holds. Hence, m and 
X scale with L as before, albeit with different exponents. Similarly, for the scaling of 
the disconnected susceptibility at Tc, we expect that 

XL.dis cx L^'\ (20) 
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If we assume that each distribution PL.i{s) scales conform Eq.®, it follows that 
{\s\)i and {s'')i scale according to Eqs. ljlOaP and (jlQ&p . respectively (the subscript 
denotes that the thermal average was taken in the i-th random field realization). Of 
course, the scaling functions / may depend on the particular random-field realization, 
but the leading L dependence will be the same each time. Since by definition 
[{X)] = {^/M)J2i=i{X)i, it follows trivially that the L-dependence implied by 
Eqs. l|10aP and (|10&p . appears in the quenched average also. We thus obtain 

with redefined scaling functions c, which can be expressed in terms of the functions /, 
of course, but for our subsequent discussion the precise form does not matter. Using 
the definitions of x, x' and Xdis the above equation implies 

X = (C2 - ci) L'^-^'^/^ x' = (52 - Co) i'^-^''/^ Xdis = aiL'^-2^/^ (22) 

On the other hand, finite size scaling also demands that x o^ x' oc L''^" and 
Xi,dis oc L^''^ . The solution of the paradox is to require that 

7 + 2/3 = i/d, (23) 

which correctly sets the scaHng of Xdis, and also that Cq = £2 and Ci = 62. Note 
that Eq. (|23|) is just the standard hyperscahng relation, but with 7 replaced by 7. 
Hence, even though normal hyperscahng in the RFIM does not hold, Eq.® still gives 
a consistent description of finite size scaling, but one must accept that the connected 
susceptibility is not described by it, since the leading terms in x and x' cancel. To 
also describe the scaling of x and x\ one needs to include the leading correction to 
scaHng. This correction can be derived by assuming that PL,i{s) at and below Tc is 
a superposition of two Gaussians. Expressing the peak at positive magnetization as 
P^iis) oc exp (-(s - 771^)^/(2^1)), it follows that (s)i — m^ and (s^}i ~ mf +w1. 
Performing the quenched average, we now obtain a non-zero expression for the 
connected susceptibility x = [JJ^ jM) X)i=i "^1- This term, consequently, is the sought- 
for correction; finite size scaHng then implies that wf ex ljil^-<^ at criticality. 

Since cq = £2, it also foHows that the cumulant at criticality U* = [(s^)]/[(|s|)^] in 
the RFIM tends to unity [42,43]. The shape of the quenched-averaged distribution at 
Tc is therefore similar to that below T^: both distributions are characterized by C/i = 1 
in the thermodynamic limit. Hence, also at Tc, we have a distribution featuring two 6- 
peakslj. This is profoundly different from systems where hyperscahng holds, since here 
Ui tends to a non-trivial value different from the off-critical values (as explained in the 
previous section). For the RFIM, plots of C/i versus T, for various system sizes L ^ 00, 
no longer intersect. In practice, however, the system sizes feasible in simulations are 
still quite small, and so one is plagued by "cross-over" effects [52] (in this case from Ising 
to RFIM universality). This means that an intersection point can typically still be 
identified, but it occurs at a value much closer to C/j^ = 1 of the RFIM [37,43,50][|J. The 
fact that the quenched-averaged distribution in the RFIM remains sharp at criticality, 
in contrast to overlapping, is also obvious from Eq. (|23|) . Considering again the ratio 
A between peak- width and peak-to-peak distance, i.e. conform Eq. (fT3|) . one finds that 
u> — (7 — 7)/(2i'). Using the result of Schwartz that 7 = 27 [45], it immediately follows 
that uj < 0. In other words, for the RFIM at its critical point, the relative peak width 
A vanishes, leading to a distribution featuring two (5-peaks. 

§ Of course, whereas for T = Tc the peak positions scale oc L"''/", they saturate at finite values ±m 

when T <Tc. 

II Note in particular Figure 7 of [50] for the RFIM. 
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Figure 2. (a) {s) versus H at T > Tc, for one realization of ttie random field 
where He happens to be positive. At H = 0, we have a finite magnetization 
of order xhL~'^/^ (point A). The field at which {s) changes sign, and where 
X ot d{s)/dH attains its maximum, occurs at -ff = —H^ (point B). Note that the 
slope at B approaches the zero-field connected susceptibility x\h=o '° ^^^ limit 
L ^ oo. (b) [(s)] versus -ff at T < Tc. The dashed curve shows the behavior 
in the thermodynamic limit; the solid curve in a finite system of size L. Note 
that the rounding is of order L""^'^. This means that the slope d[{s)]/dH\g^Q in 
finite systems grows oc L'*'^, and that the region where [(s)] deviates significantly 
from L — > oo behavior shrinks oc L~ ' . 



2.3. random-field Ising model: sample-to-sample fluctuations 

The result of Schwartz [45], namely that 7 = 27, can be made plausible when we 
consider one particular reaHzation of the random field. In a volume L'', roughly half 
the lattice sites "feel" a negative random field (and the other half a positive random 
field, obviously) but with Poissonian fluctuations. Hence, there will typically be an 
excess Zeeman energy of order ±/iL'*'^, which has the same physical effect as if an 
uniform external field of strength 

He - ±hL-'^/^, (24) 

acted on the spins in this volume (recall that h is the strength of the random field). 
But then we expect a non-zero magnetization (s) = HcX ^ x^L~'^/'^ in this sample, 
with X the connected susceptibility. Using near Tc the standard finite size scaling 
relations for (s) and x, we obtain L~f^/^ oc LP'I^^'^I'^ ^ or fijv = d/2 — ^/v. Combining 



with Eq.|[23l) one finds that 7 = 27. 

It is of some interest to explore the consequences of Eq. (|24|) further, and study the 
behavior of the magnetization for different realizations of the random field. If TJc > 
and T > Tc, we have for if = a positive magnetization of order (s) ^ x/iL~'*/^ as 
argued above (recall that H is the strength of the uniform external field) . The field at 
which the susceptibility x '^ d{s)/dH is maximized is therefore not H = 0, but rather 
H = —He, where the net effect of the random field is canceled. This is sketched in 
Fig. [2l[a), where (s) versus H is plotted [31. Taking the limit L — > co, it follows that 
the slope of the curve at points A and B becomes the same, since, on the small scale of 
He, the curve may be approximated by a straight line. Note that the slope approaches 

^ Of course, the graph of [(s)] versus H is anti-symmetric about the origin, since in the quenched 
average both signs of He appear equally often. 
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(a) pure Ising model 
T<T 




*-H 



Figure 3. (s) versus H aXT <Tc for the pure Ising model (a) and the RFIM (b). 
The dashed curves show the behavior in the thermodynamic limit; solid curves 
for finite systems (see details in text). Note that scenario (b) holds only in spatial 
dimension d > 2. 



the zero-field susceptibility x, and also that the slope is independent of He. Plotting 
(s) versus H for different realizations of the random field, one thus obtains a set of 
parallel straight lines. In other words, x is rather insensitive to the particular random 
field configuration, which just expresses the fact that the system is self-averaging for 
T > Tc, as expected [53,54]. This result is important since, in systems lacking spin 
reversal symmetry, the natural path in the (T, _ff )-plane to follow is no longer the fine 
H = 0, but rather the path along which x assumes its maximum for each realization 
of the random field. 

The situation is qualitatively different for T < Tc, of course, since we now 
expect a first-order transition and, consequently, two-phase coexistence. In the pure 
Ising model, coexistence between two states (with positive and negative spontaneous 
magnetization ±m) occurs at H = 0, irrespective of the system size L. In the 
thermodynamic limit, (s) increases monotonically with H, and jumps from —m to 
+m at if = 0; see the dotted curve in Fig. Elja). In a finite system, the transition 
is rounded, and a true jump does not appear. Instead, (s) passes smoothly through 
the origin, but with slope d{s)/dH oc L'^ [55,56], consistent with the formation of a 
jump in the thermodynamic Hmit; see the full curve in Fig.jSja). Consequently, phase 
coexistence in the pure Ising model may always be studied using H — 0. Provided 
T is sufficiently below Tc, double-peaked distributions Phis) are readily observed, 
i.e. conform Fig. [TJa), from which the coexistence properties follow. 

For the RFIM below Tc and finite system size L, the behavior is more subtle, 
since the random field breaks the spin reversal symmetry. We still expect a (rounded) 
first-order transition, but centered around the shifted field H — —He, see Fig. [3ljb). 
In the thermodynamic limit. He — > 0, and so the magnetization jumps, as before, 
at H = (dotted curve). In a finite system (s) increases smoothly with H (solid 
curve), passing through zero at H = —He (point A), with slope d{s)/dH oc L''. Since 
He ex hL^'^/'^ asymptotically exceeds the rounding, it follows that, in a finite system 



Colloid-polymer mixtures in random porous media 11 

at H ~ 0, phase coexistence is unlikely. At H ^ 0, one either observes the phase with 
positive magnetization (as one would in Fig. [H^b)), or, if the random field happens 
to resemble He < 0, a negative magnetization. Only very rarely, when the inflection 
point A happens to coincide with H = 0, will both phases be observed simultaneously. 
Hence, at H — 0, the distribution PlaIs) will mostly feature just one peak, located 
at positive or negative values. In the quenched average, one recovers [(s)] = 0, of 
course, but the fluctuation [(s)^] — [(s)]^ is not zero, since this, apart from a factor 
L"^, is precisely the disconnected susceptibility, see Eq. (|15cp . Clearly, to study phase 
coexistence in simulations, it does not make sense to use H — 0, since one would 
rarely see a double-peaked distribution. Instead, it is more meaningful to obtain these 
properties at the inflection point A, where x oc d{s)/dH attains its maximum. To be 
precise: one should apply an external field H — —He "tailored" for each random-field 
reaHzation. In the Hmit L — > oo, one has He -^ 0, and coexistence properties obtained 
at the inflection point, will agree with those obtained at H — 0. The advantage of 
the former method being that double-peaked distributions Pla{s) will now already 
appear in much smaller systems. Of course, for these double-peaked distributions, 
{s)i will be close to zero each time, and so it follows that a different definition for the 
disconnected susceptibility should be used, presumably of the form Xj;,, = L'^[(|s|)^]. 

Considering now the behavior of [(s)] versus H below T^, we expect the scenario 
of Fig. [2l[b). In the thermodynamic limit, a jump in [(s)] at iJ = is anticipated. In 
finite systems, the jump is rounded, but on a more severe scale i"'''^, as pointed out 
by Kierlik et al [14]. Note that graphs of [(s)] versus H for finite L intersect the origin 
since, in the quenched average, both signs of He are equally likely. 

Unfortunately, these arguments cannot be easily extended to T = Tc. As 
discussed in detail by Wiseman and Domany [53,54], systems with quenched random 
disorder at criticality exhibit lack of self-averaging. For small enough fields H , it still 
holds that (s) versus H for one realization of the random field, is a straight line, with 
slope oc L''/" . The same holds for [(s)] versus H, where the slope is also oc L'^^" , but 
the prefactors differ. The ratio of these prefactors is a quantity characterizing the lack 
of self-averaging, in the sense of Wiseman and Domany [53,54]. 

2.4- obtaining the quenched average using a sample dependent Hi 

For the Ising model, one knows beforehand that the infiection point of (s) versus H 
(or [(s)] versus H in case of the RFIM), occurs on the symmetry line H — 0. Hence, 
varying T at fixed H — 0, one cannot miss the critical point. In less symmetric 
models, the field H at the infiection point is not known beforehand. In these cases, 
it is clearly more convenient to follow the path H — —He{T) in the (T, iJ)-plane 
of each random-field realization. That is, for each realization of the random field i, 
one numerically locates the field Hi where d{s)i/dH in that sample is maximized. 
Properties of interest are then collected at Hi , and the process is repeated over many 
different random field samples. Extrapolating the data to L — > oo is demanding in 
practice, but does not present any principal objections. Only the prefactors of the 
finite size scaHng laws at criticality 

^cxL"^/^ (s2)-(s)2 o^L-'/"-'^, l/p o^L^I"-'^, (25) 

will differ from those of the standard quenched-average [•] obtained at fixed H. Here, 
the overbar denotes averaging at the sample dependent Hi . Since many typical fiuids 
(including the AO model) are asymmetric, collecting the quenched average as X simply 
becomes a necessity in these cases. 
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Figure 4. Schematic representation of ttie expected bimodal form of PL,i{N) 
obtained in a single realization of the porous medium at and below Tc. The 
average of the full distribution (A'')^ separates the vapor from the liquid peak. 
The distance between the peaks, defined as the average number of particles in the 
liquid phase (Af)iiq_i minus the average (Af)vap,i of the vapor, gives an estimate 
of the order parameter. The moments of the vapor and liquid peak yield the 
connected and disconnected susceptibilities, see details in text. 



2.5. extension to fluids 

We now consider a liquid-vapor transition of a fluid conflned to a quenched porous 
medium. We use the grand canonical (GC) ensemble, i.e volume L'', temperature T, 
and chemical potential jjl are flxed, but the number of particles N in the system 
fluctuates. Our analysis is based on the (normalized) distribution 

PL,^{N)^PL.^{N\T,^ii), (i = l,...,M), (26) 

deflned as the probability to observe a system containing N particles, in the «-th 
realization of the porous medium. Note the dependence on L and T, and also that 
we allow the chemical potential /i^ to vary between different realizations of the porous 
medium. For given L and T, PL.i{N) is sampled from A^ = to A'max, using a 
biased sampling scheme [57]. This process is repeated for M different realizations 
of the porous medium. The sampling scheme is constructed to visit the full range 
< A^ < Njnax irrespective of the imposed chemical potential. Hence, we set ^^ = in 
the simulations, and use histogram reweighting [58] to extrapolate to different values 
afterward. 

Assuming that the liquid-vapor transition in a porous medium belongs to the 
universality class of the RFIM, we expect, in the thermodynamic limit, a critical 
point at temperature Tc and chemical potential ficr- Below Tc, we anticipate bimodal 
distributions PL,i{N), but only if /Xj is chosen reasonably. In flnite systems, we actually 
expect the bimodal form to persist considerably above Tc also, since, for the RFIM, 
PL,i{N) remains sharp at criticality. In this work, /i^ is tuned for each realization 
of the porous medium such that d{N)i/d^i for that realization is maximized, with 
{N)i = Y^^={)^PL,i{N)\^- Loosely identifying {N)i ^^ (s), ^i ^^ H, our choice may 
be regarded as the analogue of the inflection point A in Fig. [3ljb) for the magnetic 
case. 

+ Note that other choices are conceivable also, such as the "equal-area-rule" [56,59], or the generalized 
fc-locus defined in [60]; all become identical in the limit L -^ oo, of course. 
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Fig. [4] shows a schematic of PL,i{N) in bimodal form. The left peak represents the 
vapor, the right peak the Uquid, with the average {N)i of the full distribution located 
in between (arrow A) . Fig. |4] also shows that A^max should be chosen well beyond the 
liquid peak (arrow B). If we shift PL,i{N) by its average, we approximately recover 
the Ising symmetry PLis) = Pl{—s) of the magnetization distribution. Taking the 
quenched average, this requires a shift over [{N)] — (1/M) J2i=i{-^)i- Therefore, 

iN-[{N)])/L^, (27) 

in a fluid with quenched disorder, is the analogue of s in a magnetic system, where 
the factor L'^ is needed because s is the magnetization per spin. Replacing s in the 
definitions of x and Xdis by Eq. (|27|) . one obtains 

X = [{N') - {Nf]/L^ (28) 

Xdis=([(iV)2]-KA^)nM (29) 

as the analogues of the connected and disconnected susceptibility in a fluid with 
quenched disorder. As stated before, x ^ind Xdis are analyzed for the vapor and 
liquid phase separately, using {N)i as a "cut-off" separating the peaks in PL,i{N). In 
this way, we obtain for the vapor phase 

W. 
(A^^)vap,» = 2 J2 N'^PlAN), (30) 

where the factor-of-two is a consequence of the normalization of PL.i{N). The 
moments (Af'^)iiq,i of the liquid are obtained similarly, with the summation from 
N = {N)i to A^niax- The connected and disconnected susceptibilities of the vapor 
phase can now be written as 

,,, _ [(A^^)vap] - [(A^);ap] vap _ [(^)^ap] " [JN)..,? 

Xcon Td ' '^dis Td ' ^ -^ 



k\l 



vap J 



with the quenched average [•] conveniently expressed in terms of Eq. lpO]) as [(A^*"') 

(1/M) J2i=i (-^'^)vap,i- Similar expressions hold for Xcon ^^d x^? also. Note that, since 
the chemical potential fii is "flne-tuned" for each realization of the porous medium, the 
quenched average obtained above actually corresponds to X of Eq. l(25ll , but this should 
be obvious from our discussion. For completeness, we remark that the quenched- 
averaged distance between the peaks in Fig. |4] may be used as order parameter 
m = [(A^)iiq] — [(A'')vap], although in this work the emphasis is on the susceptibilities. 
Of course, it needs to be verifled in simulations if the expected bimodal form 
of PL,i{N) really occurs in practice. In our previous work, this turned out to be the 
case [37]. However, GC simulations of the Lennard- Jones fluid with quenched disorder 
have revealed distributions with three peaks also [9]; the possibility of two fluid phase 
transitions occurring has also been suggested [10], although this probably does not 
survive in the quenched average [11]. 

3. Model and simulation method 

3.1. AO model 

We now proceed to test the concepts of the previous section in a colloid-polymer 
mixture with quenched disorder. Our primary aim is to measure the connected 
and disconnected susceptibilities, and to show that both diverge at criticality. To 
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describe the mixture, we use the AO model [20,21]. In this model, colloids (species c) 
and polymers (species p) are treated as spheres with respective diameters Cc and 
Up. Hard sphere interactions are assumed between colloid-colloid and colloid-polymer 
pairs, while the polymer-polymer interaction is taken to be ideal. In this work, CTc is 
the unit of length, the colloid-to-polymer size ratio q = a^ja^ — 1, and the spatial 
dimension will be d = 3. The behavior of this model for g = 1 without quenched 
disorder has been studied before [36], and bulk phase separation, whereby the mixture 
"splits" into a colloid-rich (polymer poor) and colloid-poor (polymer rich) domain, 
was readily observed. If one "identifies" the colloid-rich phase with a liquid, and the 
colloid-poor phase with a vapor, the phase separation can be treated in much the same 
way as a liquid-vapor transition. In the GC ensemble, one then introduces the colloid 
chemical potential ^, and, following convention, the polymer "chemical potential" jypLj. 
Phase separation occurs at the coexistence colloid chemical potential ^ = /Xcoex, for 
values of ryjj exceeding the critical value rj^ ^.^ (77p is therefore the analogue of inverse 
temperature; for q = 1, r]^ ^^ fv 0.861 has been reported [36]). The discussion and 
definitions of Section 12.51 thus trivially "carry-over" to the AO model if one identifies 
A^ ^^ number of colloids, /i ^^ colloid chemical potential, and T ^^ l/'7p- 

3.2. AO model with quenched disorder 

To study the AO model with quenched disorder, we introduce a third species Q 
of immobile (quenched) particles. These particles are also spheres, with diameter 
cq = Cc , and they are distributed in the simulation box at the start of each simulation 
(the simulation box, incidentally, is a cube of volume V = L'^ with periodic boundary 
conditions). The quenched particles, A^q of them in total, are located at random 
positions, irrespective of overlap. Consequently, the structure of the quenched system 
is just that of an ideal gas. The average packing fraction of the quenched system 
is fixed at rjQ — TTaQNQ/{6V) = 0.05, but, consistent with our GC approach, we 
allow for Poissonian fluctuations around the average. Prom a computational point of 
view, the quenched system is trivial to generate: one simply draws A^q from a Poisson 
distribution, and generates a corresponding number of positions in the simulation box. 
Next, a GC simulation of the AO model is performed in the simulation box containing 
the quenched system, whereby the colloid and polymer positions are continuously 
updated, but not the positions of the quenched particles, of course. The colloids and 
polymers interact with the quenched particles in a simple way: colloids may not overlap 
with quenched particles, while the polymers may overlap freely with them. Of course, 
computational efficiency is the main motivation for using such simple interactions, 
although one could envision similar interactions in experiments also, using polymer 
quenched disorder. In any case, the simple approach adopted here is appealing, as 
previous work indicates [13,37,38]. An estimate rj^^j. « 1.192 has also already been 
reported [37], for the exact same parameters as considered here. 

3.3. implementation details 

We now discuss some implementation details. For the i-th realization of the quenched 
system, grand canonical MC is used to measure PL,i{N) of Eq. l(26l) . with A^ the 

* Strictly speaking, ri^ is defined as the polymer reservoir packing fraction [22]. For the present 

case of ideal polymers r/p = nape^'^p''''^'^' /6A^ , with /ip the polymer chemical potential, and A the 
thermal wavelength. 
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number of colloids. The distribution is obtained using the (already mentioned) biased 
sampHng scheme [57], in conjunction with a cluster move [25,61]. The cluster move is 
needed to alleviate the otherwise (too) slow equihbration of the AO model. Of course, 
simulations of a single-component fluid do not require the cluster move. We consider 
system sizes L = 7 — 12. For each system size, P^^i^N) is typically measured for 
M ~ 2000 (!) reaHzations of quenched disorder, at several values of rj^ in the vicinity 
of r7p_ci- Large values of M are needed to obtain Xdis accurately. 

In GC simulations, particles are continuously inserted and deleted from the 
simulation box, and so one can deflne a time r after which a given population 
of particles has been completely "updated" by new ones. The duration of a GC 
simulation may therefore be expressed in units of t. In the biased sampling scheme [57], 
simulation time can be conveniently allocated, since the scheme constructs PL,i{N) 
step-by-step via so-called windows. In the first window, N varies between and 1, in 
the next window between 1 and 2, and so forth, up to A^max (the number of polymers 
A^p fluctuates freely in each window, of course). Hence, we allocate a flxed amount of 
simulation time, typically 5t, to each window. It then takes roughly 12 minutes to 
obtain P^^iiN) for L = 7, and about 1 hr for L = 12. Of course, these benchmarks 
depend on rf^, as well as on the precise computer architecture, but they suffice to give 
an overall impression of how much computer time was used. 

A final remark concerns the implementation of histogram extrapolation [58]. As 
stated earlier, all simulations are performed at colloid chemical potential /i^ = 0, 
and P{N\^i — fi') oc P{N\fii = 0)exp{fi'N) is used to extrapolate to different 
values. Obviously, a similar expression holds for the polymers also, which one could 
use to extrapolate in rjp. In fact, an important ingredient of this work is precisely 
the latter extrapolation, and our analysis would become extremely cumbersome 
without it. However, this requires that we store the full two-dimensional histogram 
PL,i{N, Np), with TV the number of colloids, and A^p the number of polymers. Since we 
typically consider 2000 realizations of quenched disorder, storage requirements become 
enormous. Fortunately, storage can be reduced tremendously, when one realizes that, 
for a fixed number of colloids N, the corresponding distribution in A'p is to a good 
approximation a single Gaussian peak. For N = this is obvious, since then we 
have a pure polymer system, but it holds well for A^ > also. Hence, to facilitate 
extrapolations in r]^, we only need to store the average and variance in A^p for each 
window (which costs only very little storage, at no cost in CPU time either). We have 
verified this approach and checked that results obtained at one value of r]p indeed 
extrapolate to those obtained at a different value (not too far away, of course) . Note 
also that the histogram extrapolation method itself can be optimized since, for a 
Gaussian distribution, integrations over A^p can be performed exactly beforehand; the 
resulting expressions become functions of the average and variance, which can be 
hard-coded. 

4. Results 

4-1. sample to sample fluctuations 

The analysis of Section 12.51 requires that the distributions PL,i{N) are somewhat 
bimodal, i.e. that they resemble the schematic shape of Fig. [H In order to verify 
this, we show, in Fig.IU PL,i{N) for a number of realizations of quenched disorder, at 
a value of rj^ significantly below the critical value rj^^j. « 1.192. Clearly, the bimodal 
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Figure 5. Distributions PL,i{N) for 18 different realizations of quenched disorder 
using rjp = 1.05 and L = 10. The horizontal axes in each of the plots show the 
colloid packing fraction rjc = iTa^N/{6V) from r/c = — » 0.2 (left to right); the 
unit on the vertical axes is arbitrary. 






Figure 6. Same as Fig. [5] but for rf^ = 1.15 and L 
fraction on the horizontal axes is from r;c = — > 0.25. 



10; the colloid packing 



Colloid-polymer mixtures in random porous media 



17 




0.75 



500 



1000 



1500 



2000 





2000 



0.40 



0.40 
0.30 
0.20 
0.10 



500 



1000 



1500 



2000 



XZ for L=7 and Ti;,= 1.00 



500 



1000 



1500 



2000 



Figure 7. "Moving average" of ttie susceptibilities Xcon (left frames) and Xjii 
(right frames) of the liquid phase, for L = 7 and 11, using several values of »7p as 
indicated in the labels. Plotted are the susceptibilities (vertical axes) versus the 
number of quenched disorder realizations M (horizontal axes). 



shape is already present in most distributions, even for this low value of rj^. Of course, 
by making rj^ even lower, the bimodal shape will eventually vanish for all realizations 
of quenched disorder, since we then enter the one-phase region where PL,i{N) is just 
a single peak, conform Fig. [ijb). In any case, Fig. [5] does confirm our expectation 
that, for random-field Ising universality, bimodal distributions persist well above Tc 
(recall that rj^ is the analogue of inverse temperature). Fig. [5] also reveals that not 
all the distributions are bimodal, see, for example, the distribution in the upper left 
corner. In these cases, splitting the distribution in half at the average is not meaningful 
anymore, although numerically this can still be applied. Since, for the thousands of 
distributions generated in our simulations, inspecting each one visually by hand is not 
feasible, the (occasional) single-peaked distribution is treated in the same way as the 
bimodal ones. Of course, single-peaked distributions become increasingly rare upon 
increasing rj^, as Fig. [6] clearly indicates, where the same realizations of quenched 
disorder were used as in Fig. [H Note that, in Fig. [6l all distributions now feature two 
peaks. 

Another feature that emerges from these figures is that the vapor peak is much 
sharper than the liquid peak. This appears to be a non-universal feature that depends 
on the interaction between fluid and quenched particles. In our previous work, we 
have studied a different type of quenched disorder, whereby also the polymers were 
not allowed to overlap with the quenched species [37,38]. In this case, a reversed trend 
was observed, namely a sharp liquid peak "coexisting" with a much broader vapor. 
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Figure 8. (a) Connected susceptibility of ttie liquid phase versus rj^ for several 
system sizes L. Note the increase of peak height with L, and also the shift in 
the peak positions, (b) Finite size scaling plot, where r]^^^^. = 1.194, u = 1.1, and 
■y/u = 1.87 were used (see details in text). 



Having shown that our assumption of bimodal distribution shape is a reasonable 
one, we need to determine the number of quenched disorder reaHzations M typically 
required to obtain Xcon ^^^ Xdi^ accurately. To this end we show, in Fig. [7l the 
variation of these quantities as a function of M, for two system sizes L, and several 
values of rjp (as indicated in the label of each subplot). The trends revealed in Fig. [7] 
are typical for other state-points also. Clearly, from this figure, we conclude that M 
should be of the order of 1000 at least. Larger values are better still, but then we 
meet the limit of our computational resources. 



4-2. connected susceptibility 

We now consider the connected susceptibility, first of the Hquid phase. Shown in 
Fig. EJa) is Xcon versus rj^, for several system sizes. Note the presence of the peak. 
Consistent with finite size scaling, the peak height increases with L; the latter could 
now be fitted to Eq. lfffcl) to obtain j/i^. However, a more stringent test is to plot 

tL'/'^ vs. xtconL--^'"^ (32) 

with t = ?7p/?7p.cr ~ 1 the relative distance from the critical point. Although not 
derived in this work, finite size scaHng implies that data from different system sizes, 
when scaled conform Eq. (|32|) . collapse onto a single master curve, provided the correct 
values of rf^ctJ ^i ^^d 7 are usedljj. In Fig. [8l[b) the resulting scaHng plot is shown, 
where rf ^^ — 1.194, v — 1.1, and 7/1^ = 1.87 were used. The quality of the collapse is 
clearly very good. However, we noticed that good collapses were obtained for different 
values also, typically i/ = 1.0 — 1.2 and r7p ^.j. = 1.19 — 1.22, which gives an indication of 



The derivation is straightforward, see for example [41]. 
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Figure 9. Same as Fig. [8] but for ttie connected susceptibility of t tie vapor. In 
tiie scaling plot of (b), »?p_cr = 1.194, ly = 1.1, and ■y/u = 1.87 were used. 



the uncertainty. The problem is that both v and rjp^^^ follow from the X-dependence 
of the peak positions. Over the range of available system sizes, the shift in the peak 
positions is rather small, and hence large uncertainties in v and 77p ^.^ are unavoidable. 
In contrast, 7/1^ can be obtained more reliably, since the latter is set by the peak height 
versus L, which yields a more pronounced numerical signature. Similar conclusions 
are reached for the connected susceptibility of the vapor, see Fig. [9l 



4.3. disconnected susceptibility 

We now come to the main result of this work, namely the behavior of the disconnected 
susceptibility. If fluids with quenched disorder belong to the universality class of 
the RFIM, the analogue of Xdis defined in Section 12.51 should diverge with critical 
exponent 7. Since 7 = 27 [45], the divergence should be very pronounced, much 
more pronounced than that of the connected susceptibility, in fact. In Fig. [TOF a). we 
show Xdii of the Hquid phase versus rj^ for several system sizes. The formation of a 
peak is clearly visible. Note also the rapid growth of the peak height: increasing the 
system size from L — 7 ^ 12, the disconnected susceptibility increases by a factor of 
more than six, compared to a factor of about three for the connected susceptibility. 
The corresponding scaling plot is shown in Fig. [TOlfb). which now involves 7, of 
course. Using rj^,,^ = 1.194, i/ = 1.1, and j/v = 3.82, the data collapse convincingly, 
confirming the power law divergence of Xdis- For the same reason as before, the scaHng 
plot is rather insensitive to ryp ^.^ and v, and so the uncertainty in these quantities is 
similar as before, but the ratio 7/1^ should again prove reliable. In Fig. HDJa) we plot 
the disconnected susceptibility of the vapor phase versus rj^, but only for L < 10. 
For reasons we do not yet fully understand, the statistical uncertainty in Xdts'' ^^ ^^'^y 
large. While the growth of a peak with system size is still confirmed, the data clearly 
do not lend themselves for measuring critical exponents, and so a scaHng plot is not 
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Figure 10. The main result of ttiis paper: (a) disconnected susceptibility of the 
liquid phase versus rj^ for several system sizes L, and (b) the corresponding finite 
size scaling plot, where »?p,cr = 1-194, u = 1.1, and ^/u = 3.82 were used (see 
details in text). 
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Figure 11. (a) Disconnected susceptibility of the vapor phase versus Jjp, and (b) 
the quenched-averaged chemical potential versus r^jj (for several system sizes L). 
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Figure 12. Plots of Ui versus ri^, using several system sizes L, for the AO model 
with quenched disorder (a), and without (b). The horizontal line in (b) marks 
C/j ~ 1.239 of the Ising model in three dimensions. 



attempted here. One reason for the large statistical uncertainty in Xdis i^ the smaller 
number of particles in the vapor phase (compared to the Hquid) . 



4-4- scaling of the cumulant 



In Fig. [T2i a). we show the cumulant as a function of r^p for several system sizes. Recall 
that the cumulant is defined as Ui = [(s^)]/[(|s|)^], with s ^ {N - [(N)]) /L'^, which 
can be calculated straightforwardly from the distributions PL,i{N) U. As expected, 
the cumulants from different system sizes do not intersect at criticality, but instead 
reveal a scatter of intersection points, close to U* — 1 of the RFIM. This behavior is 
conform our discussion of Section [521 and confirms that PL,i{N) remains sharp at the 
critical point, featuring two well-separated peaks, since hyperscaling is now violated. 
For comparison. Fig. [T2lfb) shows the cumulant of the AO model in the pure 
system, i.e. without quenched disorder. In this case, hyperscaling is not violated, and 
a sharp intersection point is indeed revealed, occurring at a value Ui different from 
the off-critical values 1 and 7r/2, respectively. For Ising systems in d = 3 dimensions, 
we expect that U^ « 1.239 [62], marked by the horizontal Hue in Fig.[T2lfb). and our 
data indeed intersect close to this value (some deviation is clearly apparent, but to 
account for this would require a field-mixing analysis [26,63]). From the intersection 
point, we also conclude rf^^r ~ 0.876 for the pure system, which compares well to the 
estimate reported in [36]. 



I In our previous work [37], we used U\ = [(s^)]/[(|s|)]^, 
that both definitions become equivalent for L ^ oo. 



but the reader can verify following Section [2T2l 
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Figure 13. (a) Variance [/.t^] — [fi]^ versus r/^ for system sizes L = 7 — ► 12 (from 
top to bottom), (b) Same as above, but witii tiie variance scaled by L'^. 



4-5. chemical potential 

Finally, we consider the average and variance of the chemical potentials at which 
our data were obtained. Recall that, for each realization i of quenched disorder, we 
use a "fine-tuned" chemical potential fii, chosen at the maximum of d{N)i/dfj.i for 
that realization. Hence, it is interesting to consider the quenched-averaged chemical 
potential [fj] and its variance [fM^] — [/i]^, with [/i*^] — {^/M)J2i=ilJ-i- Shown in 
Fig. fTlT b) is [/J,] versus rj^ for several system sizes. The data do not reveal any strong 
L-dependence, which is similar to that observed in fluids without quenched disorder. 
Of more interest is the variance, which should vanish for L — > cx). Shown in Fig.fTSlfa) 
is [fi^] — [/i]^ versus r]^, and the decrease of this quantity with increasing L is clearly 
visible. Kierlik et al have shown that, below Tc in the two-phase region, the variance of 
the chemical potential vanishes oc L"'' [14]. Plotting therefore L'^ ([/i^] — [^]^) versus 
rip, see Fig.fTsFb). we observe that this prediction holds quite well for our data also. 



5. Discussion and summary 

We have explained finite size scaling in the random-field Ising model, and shown 
how this technique may be applied to a fiuid with quenched disorder. We have 
also defined the analogue of the disconnected susceptibility Xdis for the latter. If 
fiuids with quenched disorder belong to the universality class of the random-field Ising 
model, as conjectured by de Gennes [7], Xdis should diverge at criticality, and so our 
definition facilitates further tests of this conjecture. To perform one such test has been 
the topic of the present work, using the Asakura-Oosawa model of a colloid-polymer 
mixture confined to a random porous medium. Our data are indeed compatible with 
a divergence of Xdis- Moreover, for the liquid phase, we even recover 7 « 27, in 
quantitative agreement with the prediction of Schwartz for the random-field Ising 
model [45]. Our estimate of the correlation length exponent ly « 1.0 — 1.2, although 
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not very precise, is also consistent with reported random-field Ising estimates [50,51]. 
Hence, the present results confirm our earlier study [37], where evidence of random- 
field Ising universality in fluids with quenched disorder was also presented, but based 
on the shape of PL,i(,N) at criticality. We also remind the reader of the large number 
of porous medium realizations used in our analysis. As the "moving averages" of Fig. [7] 
indicate, such numbers become a necessity, if % and Xdis are to be obtained with any 
meaningful accuracy. 

Finally, we turn to a discussion of possible appHcations of our work to experiments. 
The prototype experimental realization of a fluid with quenched disorder is an atomic 
fluid injected into silica aerogel. This realization has the disadvantage that the 
coupling between the porous medium and the fluid is weak [3,4], as manifested by 
the small shift of the critical temperature (compared to the system without quenched 
disorder). Moreover, the characteristic length over which the aerogel structure 
appears random is very large, compared to the size of the fluid molecules. In this 
respect, colloidal fluids may offer an attractive alternative. Note that investigations 
of critical phenomena in colloid-polymer mixtures without quenched disorder [64, 65] 
are already experimentally feasible: critical interface and density fluctuations can be 
visuaHzed directly [66,67] using confocal microscopy [68]. In principle, such confocal 
experiments could be extended to include quenched disorder also. The generation and 
synthetization of quenched colloidal porous media has received considerable attention 
[69-71]. One could envision an experiment whereby a colloid-polymer mixture is 
injected into a rigid colloidal gel. Such gels could be formed using small nanoparticles 
which can grow into randomly branched networks at volume fractions of only a few 
percent [71]. The size of these nanoparticles can be much smaller than the typical 
colloid or polymer diameter, and so one can easily reach the regime where the critical 
correlations of the colloid-polymer mixture average over the random structure of 
the gel. Another feasible realization would be to use a polymer blend containing 
nanoparticles of suitable size, such that the diffusion of these particles in the blend 
is small. The structure formed by the nanoparticles will then appear to be frozen 
(quenched) on the timescales needed for the critical correlations of the polymer blend 
to equilibrate. The latter could then be measured using, for example, light scattering!^. 
In any case, we hope that the simulational efforts of the present work will stimulate 
experimental efforts also, in order to completely settle this longstanding problem. 
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